function neglnf = mlmvnorm(D,initparam)

    %Data
    x = D(:,1); %ln DC wealth with 0's when DC wealth = 0
    y = D(:,2); %ln non-DC wealth (>0)
    wts = D(:,4); %sample weights
    
    %parameters
    mu1 = initparam(1);
    mu2 = initparam(2);
    sigma1 = initparam(3);
    sigma2 = initparam(4);
    sigma12 = initparam(5);
    
    %Compute log likelihood
    S = [sigma1^2 sigma12; sigma12 sigma2^2];
    aj = (-log(2*pi)-.5*log(det(S)));
    pj = zeros(length(x),1);
    for j = 1:length(x)
        pj(j) = -.5.*[x(j)-mu1 y(j)-mu2]/S*[(x(j)-mu1)';(y(j)-mu2)'];
    end
    lnfj = wts.*(aj+pj);
    lnf = sum(lnfj);
    neglnf = -1*lnf;
    
end